Voltage sensing in ion channels: 
Mesoscale simulations of biological devices 
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Electrical signaling via voltage-gated ion channels depends upon the function of a voltage sensor 
(VS), identified with the S1-S4 domain in voltage-gated K^ channels. Here we investigate some 
energetic aspects of the sliding-helix model of the VS using simulations based on VS charges, linear 
dielectrics and whole-body motion. Model electrostatics in voltage-clamped boundary conditions 
are solved using a boundary element method. The statistical mechanical consequences of the elec- 
trostatic configurational energy are computed to gain insight into the sliding-helix mechanism and 
to predict experimentally measured ensemble properties such as gating charge displaced by an ap- 
plied voltage. Those consequences and ensemble properties are investigated for two alternate S4 
configurations, a- and 3io-helical. Both forms of VS are found to have an inherent electrostatic 
stability. Maximal charge displacement is limited by geometry, specifically the range of movement 
where S4 charges and countercharges overlap in the region of weak dielectric. Charge displacement 
responds more steeply to voltage in the a-helical than in the 3io-helical sensor. This difference is 
due to differences on the order of 0.1 eV in the landscapes of electrostatic energy. As a step toward 
integrating these VS models into a full-channel model, we include a hypothetical external load in 
the Hamiltonian of the system and analyze the energetic input-output relation of the VS. 

PACS numbers: 87.16.A- 
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I. INTRODUCTION 

Electrical excitability of cells is possible because the 
movement of a few charges can control the flow of many 
charges. This principle — amplification — led Hodgkin 
and Huxley [l| to their theory of the action potential in 
terms of electrically controlled membrane conductances. 
Electrically controlled conductances have been localized 
to channel proteins conducting Na+, K+ or Ca^"*" ions 
across the cell membrane. Voltage-controlled ionic con- 
ductances and the controlling intrinsic charge movements 
('gating currents') of ion channels have been studied ex- 
perimentally (reviewed in '4). Electrophysiology has been 
complemented by techniques measuring channel topol- 
ogy, channel structure and the change in channel func- 
tion (reviewed iny-Q)- Together, these perspectives pro- 
vide detailed information on the 'voltage sensor' (VS) 
common to these channels and exemplified by the SI- 
S4 transmembrane segments of Shaker-type K chan- 
nels. Here we use simulation to determine consequences 
of voltage sensor models that are based on the 'sliding- 
helix' hypothesis in which the charged S4 segment re- 
sponds to changes of the membrane electrical field by 
sliding in a canal formed between other transmembrane 
domains. This hypothesis, first proposed by Catterall 
, qualitatively correlates a large body of experimental 



work (reviewed in yia) • 

The relaxation times of the voltage sensor in K"*" chan- 
nels are of the order of milliseconds and thus are highly 
averaged manifestations of atomic motion in a condensed 
hase, where momentum is scattered in a picosecond 
Furthermore, site-directed mutagenesis experiments 
have shown that only a limited number of amino-acid 
residues of a voltage-dependent ion channel are individ- 
ually important for voltage sensitivity [3, l8l-[l0{ . On the 
basis of these experimental facts, we describe the volt- 
age sensor in a mesoscale model in which atoms are not 
made explicit. Because amino-acid residues with for- 
mally charged side chains in the S4 as well as S2 and S3 
transmembrane segments strongly determine VS function 
[llj , we explicitly represent the charges of these residues 
in the model. These charges are embedded in an envi- 
ronment of spatially non-uniform electrical polarizabil- 
ity. We distinguish three dielectric regions in the model 
(membrane lipid, baths, and protein) and describe elec- 
trical polarizabilities by linear dielectric coefficients that 
are uniform within each region. 

The resulting model shares important elements with 
models previously investigated by Lecar et al. [l2| , Grabe 
|13|, and Silva et al. 
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metrical details than those models because we wish to 
start from a minimalistic model; with that approach, the 
importance of structural features can be discovered as 
features are included or varied. Additionally, we use self- 
consistent methods in solving our models. 

Electrophysiological data on voltage sensor function 



are recorded by a macroscopic setup that picks up charge 
movements of a large ensemble of ion channels while 
imposing a controlled voltage across the cell membrane 
('voltage clamp' [15|, llg]). To simulate such an experi- 
ment, we encapsulate the microscopic model of a VS with 
conductors clamped to imposed potentials from external 
sources while monitoring charge displacement at these 
electrodes. The electrostatics of the system composed 
of the VS model and the electrode setup is solved self- 
consistently. The configuration space of the model VS 
(S4 translation and rotation) is systematically sampled 
to construct a partition function based on the electro- 
static configurational energy. Using this partition func- 
tion, we compute ensemble expectations of observable 
random variables (e.g., of gating charge displaced at an 
applied voltage). In this way our theoretical results meet 
two criteria for practical usefulness: they are consistent 
solutions of the physics included in the model, and they 
directly pertain to macroscopic experimental results. 

This paper describes a building block toward creating a 
mesoscale physical model for voltage-gated ion channels. 
Such a channel comprises a gated central pore (through 
which ions flow when the pore is gated open) surrounded 
in the lipid membrane by four voltage sensor domains. 
Here we describe and test a simulation system for a sin- 
gle voltage sensor. Using this setup, we simulate an 'idle' 
voltage sensor to study its reconfiguration and energetics 
as the voltage changes while no external work is done. As 
a step toward understanding how a sensor might inter- 
act with other parts of the channel, we then test how a 
simulated 'load' alters the sensor movement and how the 
work done on the load depends on the sensed variable, the 
voltage. These simulations are done for two structures of 
voltage sensor in which the S4 helix is arranged either 
in the a conformation or the 3io-conformation. These 
helix forms are currently discussed as possible alterna- 
tives for S4 structure [17h20| |. Simulations with these 
alternate structures reveal substantial consequences for 
voltage sensing. 



II. MODEL AND BOUNDARY CONDITIONS 

Figure [T] represents the simulation cell by an axial 
cross-section of the radially symmetric three-dimensional 
domain swept by rotating that cross-section about its 
vertical axis. The external boundaries (in green, la- 
beled bath & guard electrode) are the (voltage-clamp) 
electrode surfaces kept at controlled electrical potentials. 
The blue zones represent aqueous baths (labeled with a 
dielectric coefficient e^ = 80). The pink zone is a region 
of small dielectric coefficient (labeled with em = 2) that 
represents the lipid membrane. The brown zone (labeled 
S1-S3 & S4) represents the region of the channel protein 
that we model; this region is assigned a dielectric coeffi- 
cient of ep = 4 unless otherwise noted. These dielectrics 
are piecewise uniform and therefore have sharp bound- 
aries (solid black lines). Point source charges represent- 



ing protein charges of interest are embedded in the region 
of protein dielectric. Variation of their placement is part 
of this study and will be detailed later. The protein re- 
gion as seen here represents the matrix of the S4 helix 
as a central cylinder, surrounded by the other parts of 
the channel that also create a dielectric environment dif- 
ferent from the dielectric environment of the membrane 
lipid. Included in the protein region are the invaginations 
which allow the baths to extend into the planes defining 
the lipid phase of the membrane. The radius of the S4 
dielectric domain is 1 nm (a helix) or 0.98 nm (3io helix). 

Point source charges representing protein charges are 
arranged at a minimal distance of 0.2 nm from the pro- 
tein/water boundary. The charged guanido group of each 
arginine residue of the S4 segment is represented as three 
point source charges of -|-(l/3)eo on a circle of radius 
0.122 nm [blue (dark gray) spheres in Fig. [5]. The cen- 
ters of the S4 arginine charges are arranged on a helix 
defined by arginine side chains on an a- or 3io helix back- 
bone, where every third amino acid is an arginine. For 
the a helix, charged residues are separated by 0.45 nm 
in the transmembrane direction and 60° leftward around 
the helix; for the 3io helix, charged residues are separated 
by 0.6 nm and 0°. 

Dimensions of the simulation surfaces are varied within 
those constraints as defined in Figs.[ljB) and[TJC). Fig- 
ure[UB) maps the topology of the protein and membrane 
surface of Fig. [IJA) onto metrics for simulations, defined 
in Fig. [IJC). Since the surfaces in the system are radially 
symmetrical and smooth, the system is defined by a set 
of infiection points with their curvature on the left half 
of the system to be simulated. Since the gating canal is 
symmetrical in these models, this further reduces to the 
upper half of the left side. 

We analyze two degrees of freedom for the movement of 
the S4 segment: the curve on which (triplets of) S4 charge 
centers are aligned can be both translated along the helix 
axis and rotated about that axis. The model S4 charges 
thus move like parts of a solid body. Negatively charged 
residues contributed by the S2 and S3 transmembrane 
segments of the natural VS are modeled as point source 
charges of — 1 Cq fixed on a curve parallel to the curve on 
which S4 charges are centered [red (light gray) spheres in 
Fig. [2] ; the offset from the helix axis of the countercharge 
curve is 0.466 nm larger than the radius of the curve of 
the centers of the (triplets of) S4 charges [see Fig.ITjC)]. 
The axial and angular intervals between countercharges 
are discussed in Sec. IIVI countercharges are stationary in 
their assigned positions. 

The electrodes encapsulating the simulation cell serve 
three purposes: 

1. The bath electrodes provide Dirichlet boundary 
conditions corresponding to a voltage clamp. 

2. The bath electrodes substitute for screening 
by bath ions of uncompensated protein charge. 
Screening by the ions in an aqueous bath is equiva- 
lent to the screening provided by charge on a metal 
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FIG. 1. Simulation cell. (A) The 3D setup is produced by rotating the cross-section about its vertical axis. Green lines (labeled 
as bath & guard electrodes) are electrode surfaces bounding the cell. Black lines are dielectric boundaries separating uniform 
dielectric regions: baths {blue, with dielectric coefficient Ew), membrane lipids {pink, with dielectric coefficient tm), and protein 
{brown with dielectric coefficient tp). Charges of protein side chains (represented by colored balls in Fig. [2| are embedded in 
the protein dielectric region in varied geometries. We thus simulate a single VS sensor domain (S1-S4) modeled as a central S4 
cylinder surrounded by a ring of protein material including the S1-S3 transmembrane domains. The junction between these 
protein domains is narrowed to less than the membrane thickness by circular invaginations ('vestibules') leading up to the 
'gating canal' through which the S4 helix glides through the rest of the protein (dashed line). (B, C) Geometrical parameters 
of models. (B) Mapping from geometrical positions in (A) to indexed geometrical parameters. Inflection points and lengths 
are varied among different models, depending on countercharge positions and helix conformation. Positions in r (radial) and z 
(axial) coordinates are marked by a colored point and an associated number. All corners are rounded with curvature radius of 
0.15 nm. Points 1,2,3,4, and 6 define the profile of the lipid and protein dielectrics, from the outermost end of the lipid domain 
(1) to the face of the S4 cylinder (6). Point 5 marks the radial position assigned to countercharges, point 7 the radial position 
of S4 charges. (C) Coordinates (in nanometers) of the points defining membrane and protein metrics. Models are symmetrical 
with respect to the z = plane. See Figs. [9] & [10] for three-dimensional representations of these geometries. 



foil placed in the water a distance away from the 
protein boundary. In the Debye-Hiickel theory, an 
electrode distance of ~0.8 nm corresponds to the 
physiological bath ionic strength; thus the electrode 
location shown in Fig. [TJA) corresponds to a bath 
solution in the low milliniolar range. An alternate 
configuration, a simulation cell with the bath re- 
gions omitted and the electrodes placed directly on 
the membrane and protein boundaries, would es- 
tablish screening at the Onsager limit approached 
at exceedingly large ionic strength. In this way, the 
possible range of screening effects can in principle 
be determined without including explicit bath ions 
in the simulation. 

3. At the surface where the membrane region meets 



the cell boundary, a set of guard electrodes forming 
rings around the cell maintain a graded far poten- 
tial varying, from ring to ring, between the poten- 
tials applied at the inner and outer bath electrodes. 
These guard electrodes impose at the membrane 
edge a far potential similar to that existing in a 
macroscopic system. 

Translational and rotational motions of the S4 helix 
are simulated by allowing the ensemble of S4 charges to 
slide within the dielectric domain of the protein shown in 
Fig. [1] The protein dielectric itself does not slide with the 
charges. With the protein dielectric extended far enough 
into the baths, keeping the S4 dielectric stationary has 
negligible electrostatic consequences because the dielec- 
tric of the model is uniform within the protein domain. 



Simulating S4 motion in this way reduces the compu- 
tational effort of model exploration by several orders of 
magnitude. 



III. METHODS 
A. Electrostatics 

We are concerned with the electrostatic interactions 
among charged groups of the VS protein, electrode 
charges and charges induced on sharp dielectric bound- 
aries. In solving the electrostatics, we take advantage of 
the fact that all other charges besides the point source 
charges of the protein are distributed on a few boundary 
surfaces rather than distributed throughout a volume. 
The primary task of calculating electrostatic interactions 
consists in determining the charge distributions on the 
electrode and dielectric boundaries, distributions which 
are initially unknown for a given configuration of protein 
charges and applied voltage. 



1. Computation of unknown charges 

Boda et al. 21] have described and tested a bound- 
ary element method (the induced charge calculation) 
for computing the charge distribution on the dielectric 
boundaries of a system consisting of point source charges 
and linear isotropic dielectrics with sharp boundaries. 
We include an additional electrostatic element, the 'elec- 
trode': an infinitesimally thin conductive foil charged to a 
prescribed electric potential using an external source. In 
our system, this surface is also a dielectric boundary be- 
tween the simulation cell and the dielectric surrounding 
the cell (e.g., vacuum). Thus the charge of an el ectrode 
is in part induced charge (as in the paper by Boda et alh 
and source charge (to fix the prescribed potential). In 
the following, we show how the total electrode charge can 
be calculated without needing to separately calculate its 
parts. The electrode charge calculation complements the 
induced charge calculation, making it possible to calcu- 
late all the initially unknown charges in the system. 

The spatial density of source charge p^'^^ in our system 
is composed of the point source charges q^^ of the VS pro- 
tein located at positions r^ and source charge distributed 
on the electrodes 8 with surface density (j^^'^{y £ £). 

The polarization charge density induced at any loca- 
tion r in a dielectric is 
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where e is the generally location-dependent dielectric 
coefficient, eo the permittivity of the vacuum, and E 
the electric field strength produced by all (source and 
induced) charges in the system. This relation follows 
from Poisson's equation (including polarization) and the 



constitutive relation describing polarization in a linear, 
isotropic dielectric ^211 . The first term on the right-hand 
side of Eq. dH) describes the charge induced on the sur- 
faces of the volume element at r if the element contains 
source charge. The second term describes charge induced 
in the volume element by the electric field if the dielectric 
coefficient at r has a non-zero gradient. 

For our charge calculations, we combine collocated 
source and induced charges into an effective charge for 
computing the field and potential. The effective charge 
density associated with a known source charge density 
psrc^-j.^ embedded in a dielectric [described by a locally 
uniform e(r)] is 
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The effective charge density of an electrode includes both 
contributions to the induced charge described by Eq. ([Ij , 
as well as the source charge. 

The dielectric boundaries B inside the simulation cell 
(marked in heavy hlack in Fig. (TJ all lines but the elec- 
trodes) do not carry source charge. However, the elec- 
tric field in the simulation cell induces the charge density 
CT'"'^(r) at locations r G S. This induced charge density 
is initially unknown. 

The field strength E and potential V in our system are 
produced by the superposition of the fields and potentials 
of the source and induced charges: 
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where da' is the area of the surface element at location 
r'. 

The unknown induced surface charge density on the 
dielectric boundary is related to the field strength by |2]| 



a'"^(r) 
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(5) 



where r is any location on the surface 6, Ae(r) is the 
difference in the dielectric coefficient across the dielec- 
tric boundary in the normal direction n(r), and e(r) is 
the mean of the dielectric coefficients at the boundary 
location. 





FIG. 2. VS charge positions and dielectric boundary surfaces in the a-heUcal model. Blue (dark gray) symbols: S4 charges, each 
represented as three point source charges of +(l/3)eo; Red (light gray) symbols: countercharges in the S2 and S3 segments, 
each represented as a single point source charge of — leo. The dielectric boundary surface is divided into curved tiles whose 
magnitudes are varied depending on their distance from the point source charges and local curvature. This is the surface grid 
used in solving the induced- charge calculation. The figure is drawn using perspective-enhancing features. See also Fig. [O] 



The potential at any location r on the electrode sur- 
faces £■ has a value V^'^ (r) imposed by the voltage clamp: 



V{r)^V^^{r) 



(6) 



Inserting the expression for the electric field strength 
from Eq. ^ into Eq. ([5]) and inserting the expression 
for the electric potential from Eq. Q into Eq. ([6]) yields 
two integral equations in terms of both a™'^ (r S B) and 
^effj-j. g £^y rp-j^g initially unknown charge densities on 
the dielectric and electrode boundaries are the joint so- 
lution of these two integral equations. 

To solve the integral equations, we follow the method of 
Boda et al. [21| . The surfaces B and £ are subdivided into 
curved surface elements. The unknown charge densities 
are approximated as uniform on each surface element. 
The two integral equations then become one system of 
linear equations in terms of the unknown charge densities 
of a finite number of surface elements. 

The inhomogeneity of surface charge in our simulation 
is greatest where point source charges of the VS pro- 
tein are close to a dielectric boundary or electrode. For 
computational efficiency, we vary the size of the surface 
subdivisions depending on the distance from the point 
source charges. A typical surface grid (comprising ^7000 
surface elements) is shown in Fig. [5] for the dielectric sur- 
faces in the simulation cell. The electrode surfaces (not 
included in Fig. [2|) are subdivided into ^1700 relatively 
large elements because of their distance from point source 
charges. 

The computation of the unknown charges on the sur- 
face elements involves solving a linear equation system 
in terms of as many unknowns (N) as there are sur- 
face elements. The coefficient matrix of this system 



is dense, therefore the LU-decomposition time increases 
with 0{N^). This computational disadvantage is greatly 
alleviated by the fact that L{7-decomposition of the coef- 
ficient matrix needs to be done only once for a given com- 
bination of surface geometry and dielectric coefficients. 
When the S4 charges are moved, solutions to this sys- 
tem of equations are obtained by back-substitution using 
the same LU-decomposed matrix. That 0{N'^) operation 
(back-substitution) is required for each sampled config- 
uration of VS charges but not for each applied voltage 
tested, as described later. 

Accuracy of charge calculation. The divergence the- 
orem (Gauss's theorem) states that 



(f e(r) eoE(r) • n(r) da ^ f p"''=(r) dr 
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where S is the closed surface around the volume V, da is 
the area of the surface element located at r G 5, n(r) is 
the normal unit vector for that surface element, dr is the 
volume of the space element located at r G V, and p^^'^ is 
the source charge density inside the volume. 

Gauss's theorem provides a sum-rule test for our solu- 
tion of the electrostatics that is applicable to the specific 
geometry used in a simulation. We verify the theorem 
at the dielectric boundary of the VS protein because: 1) 
that boundary is closest to the protein source charges of 
interest; and 2) the field at that boundary produces the 
largest density of induced charge in our system. The vol- 
ume integral in Eq. ([7]) yields the known algebraic sum of 
the point source charges assigned to the VS residues. The 
surface integral can be expressed in terms of the charges 
induced by the normal field at the discretized dielectric 
boundary. We consider the residue given by the differ- 
ence of the surface and volume integrals as a measure of 
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FIG. 3. Test of numerical accuracy using Gauss's theorem. 
Qcrror as defined by Eq. ((8| is plotted versus S4 translation 
[Eq. (jlSp is applied to calculate the expectation value over 
the rotational degree of fireedom] . Normal electric flux is inte- 
grated over the surface of the protein (brown region in Fig. [l] 
labeled S1-S4). The protein region contains a net charge 
of +3eo (six positive S4 charges and three negative counter- 
charges). The charge error is computed for three applied volt- 
ages ( — 100 mV, circles; mV, line; 100 mV, crosses). 



our numerical error: 



by monitoring electric flux across a surface surrounding 
a bath electrode. A more efficient method is provided by 
the Ramo-Shockley (RS) theorem [22, [23l; fo^' ^"^ appli- 
cation to ion channels see [24|). The RS theorem lays the 
groundwork as well for an efficient method of computing 
the electrostatic energy of VS configurations when the 
applied voltage is varied 125]. The RS theorem is appli- 
cable to systems containing linear dielectrics. 

The RS theorem can be formulated for the configura- 
tion of electrode potentials in this study. We apply equal 
and opposite potentials Vm/2 and — Vm/2 to the inter- 
nal and external bath electrodes to create a membrane 
voltage Vm (defined as internal minus external potential) . 
We determine the displaced charge in a simulation in two 
steps: 

1. Set all point source charges to zero and apply -1-1/2 
V at the internal and —1/2 V at the external bath 
electrode. Solve for the unknown electrode and in- 
duced boundary charges. From these charges, an 
electric potential Vb(r) can be computed for any 
geometrical location r of the simulation cell. 

2. In a simulation with the actual point source charges 
Qk present and arbitrary potentials -|-Vm/2 and 
— Vin/2 applied at the internal and external bath 
electrodes, determine the displaced charge Qk from 
the relation 



Qe 



E' 



'=p'^j 



«J - 2^ 9fc 



(8) 



Qk = -Qk 14)(rfc)/(l V) 



(9) 



where a j is the area of the protein surface element j , e j 
the dielectric coefficient on the out-facing side of aj , and 
cr"''^ the induced charge density computed to solve the 
electrostatics. 

The error in the induced charge calculation is likely to 
vary as S4 charges are moved in a simulation since the dis- 
tances between the S4 charges and the dielectric bound- 
ary vary. Figure [3] shows the error in induced charge 
calculation for the full range of S4 translational positions 
sampled in a typical simulation. The error is <0.008 eo 
of the actual net charge of 3 Cq assigned to the VS in this 
simulation. The results of our charge calculation are thus 
in good agreement with Gauss's theorem, as they must 
be. This test was performed for every simulation. 



2. Computation of charge displacement and electrostatic 
energy 

When the charges of the VS change position, the elec- 
tric flux toward one bath electrode generally increases 
by the same amount as the electric flux decreases toward 
the other bath electrode. To maintain a constant volt- 
age between the two electrodes, charge has to be moved 
externally between the electrodes. This charge is the ex- 
perimentally measured displaced gating charge. In prin- 
ciple, one can measure displaced charge in a simulation 



Note that Qk = for all geometrical positions r^ 
where Vb(rfc) = 0, and Qk varies between —qk/'2, 
and +qk/2 as the position of qk is varied from the 
internal to the external bath electrode. 

When several point source charges are in the sim- 
ulation, the total displaced charge is the algebraic 
sum of the displaced charges deflned by Eq. ^ for 
each point source charge: 



Q^Y.Q>^ 



(10) 



Step 1 is executed once for the chosen configuration of 
electrodes and dielectrics in the simulation (including di- 
electric coefficients). Step 2 is executed once for each 
configuration of point source charges. Since the displaced 
charge determined in step 2 is invariant with respect to 
the potentials applied at the electrodes, it needs to be 
calculated only once for each point-charge configuration, 
regardless of variation in the electrode potentials. 

The RS theorem also makes it possible to compute 
the electrostatic energy of a VS configuration with the 
algebraic sum of two terms: 



W = Wi+W2, 
determined by separate calculations: 



(11) 



In a simulation which includes the point source 
charges qk at positions r^, impose the potential 
Vf = on all electrodes and compute the self- 
energy 
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(12) 



where V^'^~^{Yk) is determined by Eq. (jj), exclud- 
ing self-interactions for each source charge g^.. 

2. Calculate the displaced charge Q corresponding to 
the point source charges and their positions using 
Eq. (fTU|) . For the imposed voltage V-^n calculate: 



W2 



-QV^ 



(13) 



Step 1 of this procedure is executed once for each sampled 
configuration of point source charges. Step 2 is executed 
repeatedly for each applied voltage that is tested. 

The calculations of displaced charge and electrostatic 
energy via Eqs. (TTOl) and ((TT|) are independently verified 
by computing the electrostatic energy by the path inte- 
gral of the electric force acting on the charges qu of the 
VS as those charges move from rj, to rj,: 



> 
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FIG. 4. Electrostatic energies computed via different routes. 
Energy scanned for a fixed diagonal path through the two 
dimensions (translation as shown on abscissa, rotation from 
-180° to -^180°). Results of Eq. HI]) {dots) and Eq. dH]) 
[Lines, aligned to dots at the negative end of the path). 
Applied voltage V-m'- —100 mV (dark gray), mV (black); 
+ 100 mV (light gray). 



AM^^^gfc / E(rfc)-cirfc 



(14) 



Here, the electric field is the field of all charges in the 
system except the charge qk itself, as defined by Eq. ([3|). 
Figure [H shows this control for three different applied 
voltages over a prescribed (diagonal) path through the 
translational and rotational dimensions of the range of 
S4 motion typically sampled by us. There is good agree- 
ment between the energies computed using the Ramo- 
Shockley theorem (dots) or through the path integral of 
force (lines). 



B. Statistical mechanics 

Displaced gating charge is experimentally measured 
from ensembles of channels and thus is an ensemble av- 
erage. Our electrostatic calculations yield both the dis- 
placed charge and the electrostatic part of the configu- 
rational energy for a given configuration of a simulated 
VS model. We consider whole-body movements of S4 
charge in two degrees of freedom: translation along the 
S4 axis and rotation about that axis. Our computational 
method is efficient enough to allow systematic sampling 
of this configuration space. We represent each dimension 
by 51 equally spaced grid nodes and compute the electro- 
static energy for the 2601 nodes of the two-dimensional 
space. 

The energy samples define a canonical partition func- 
tion describing the consequences of the electrostatics on 



the distribution of an ensemble in the discretized config- 
uration space: 



2-E' 



-AW.j/ksT 



(15) 



where i and j are the indices of the rotational and trans- 
lational discrete positions; AWjj is the electrostatic con- 
figurational energy of the voltage sensor at translational 
position i and rotational position j; ks is the Boltzmann 
constant; and T = 298.15 K is the absolute tempera- 
ture. The sampled rotational range is 360°, and a typical 
translational range is —1.625 nm to -1-1.625 nm relative 
to the central position of the S4 charges in a- helical mod- 
els (±2.102 nm for 3io-helical models). Restricting the 
configuration space this way is equivalent to including 
hard- wall potentials into the Hamiltonian [Eq. pT|) ]. 

The probability of a VS configuration for a particular 
translation i and rotation j is 



-\Wij/kBT 



(16) 



and the expectation (mean) value of a random variable 
X is 



{X) — y ^x^jP^j 

id 



-AWij/ksT 



(17) 



«j 
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We also determine the expectations of random variables 
over the rotational degree of freedom for a particular 
translational position i: 



{X^) = 



E,^.. 



,-AW,j/kBT 



^ g-AlVy/fcaT 



(18) 



The goal of our simulations is to study the energetics 
of the movement of the S4 segment for two degrees of 
freedom. We are concerned only with variation of the 
Hamiltonian due to changes in S4 position for this space. 
Since these simulations are done for varied settings of an 
external parameter (applied voltage), the energy at the 
reference position for AW must be invariant with respect 
to applied voltage. 

In our simulations, a voltage V^ is applied across 
the bath electrodes by imposing the potentials +Kn/2 
and — Vin/2 at the internal and external bath electrodes. 
Hence, there are locations in the simulation cell where 
the potential V due to the applied field and the displaced 
charge Qk for any point charge there are zero for any Vm 
[Eq. ©J. Since we are concerned with the position of the 
S4 segment which bears point charges qk at several loca- 
tions, we define the S4 reference position such that the 
total displaced charge there is zero. By Eq. ([T0|) and due 
to the polarity of the fields applied, there are positions 
for the S4 segment where Q is zero, and therefore the 
energy for those configurations is invariant with regard 
to Vtn. Because of the symmetries in this study, these 
reference positions coincide with the origin of the trans- 
lational axis (z = nm for any rotation (/); we choose 
= 0). 



C. Online Supplemental Materials 

Figures |1] & [10] and the associated Supplemental Ani- 
mations 1 [26| & 2 [27] illustrate VS geometry and move- 
ment for the simulations presented in this paper. They 
show the expectation(s) of position for VS charges super- 
imposed over the distribution of charge density. In the 
animation, voltage is changed in a ramp from —100 to 
+100 mV. Note that these animations are of the VS with- 
out an external load attached — as in Sec. IIV Al rather 
than Sec. HVB] 



IV. RESULTS AND DISCUSSION 

The relationship between the expectation of displaced 
charge and the applied membrane voltage in a sliding- 
helix model of an individual VS domain is shown in 
Fig. [51(A). The solid line represents the computed re- 
lation. Open circles represent the relationship experi- 
mentally observed in ShakerB K+ channels [Fig. 2A in 
[9| — the experimental charge per channel was divided 
by the number of channel monomers (4, where each in- 
cludes one VS domain). Three observations can be made 
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FIG. 5. Simulated charge-voltage relations (solid lines) com- 
pared to experiment (symbols). (A) a-helical S4 segment; (B) 
3io-helical S4 segment. The experimental relation applies 
to Shaker K"*" channels ^. Dashed lines: computed curves 
shifted to match the midpoint of the experimental curve (shift 
is 1/2 of the experimental range). 



by comparing the two relations: (1) the total amounts 
of charge that can be moved by large changes of volt- 
age are similar, ^3 elementary charges per VS domain; 
(2) the slopes of the two relations are similar; and (3) a 
shift along the voltage axis is needed to align the mid- 
point of the computed relation with the midpoint of the 
experimental relation (dashed line). 

The VS charges in the simulation for Fig. [SJA) are ar- 
ranged according to an a-helical geometry (Figs.[2]&[Hl). 
An analogous simulation in which the VS charges con- 
form to a 3io-helical geometry (Fig. rTO|) yields a differ- 
ent result [Fig. [SfB)]. With regard to the experimental 
charge-voltage relation, the 3io model yields less total 
charge movement over the tested voltage range and a 
smaller maximal slope. 

The models giving the charge- voltage relations of Fig.[S| 
involve idealized domain geometries and use the dielectric 



coefficient Cp = 4. Tfiese are initial clioices made in tfiis 
study — no variations were made to produce more real- 
istic predictions. The geometries of S4 charges were ide- 
alized to conform with two kinds of helices. The counter- 
charges (3) were arranged on a spiral (a helix) or straight 
line (3io helix) paralleling the curve of the S4 charges. 
The only parameter optimized in light of the experimen- 
tal results was the countercharge spacing (the transla- 
tional and angular intervals between countercharges). In 
the models used in this paper, countercharges are spaced 
at 2/3 the intervals of the S4 charges, so that at most one 
S4 charge lines up with a countercharge for any particular 
S4 position. 

The simulations reported in Fig.[5]indicate that voltage 
sensing by a sliding-helix is robust from an engineering 
point of view. In these generic models, either the a- or 
3io-helical S4 structure can produce volt age- dependent 
charge displacement, even though the two structures in- 
volve distinct configurations of the protein charges. The 
two forms of helix generate different responses to voltage 
in the tested models, but neither form produces catas- 
trophic failure. The biological channels whose structure- 
function relation we seek to ultimately understand in 
engineering terms themselves show robustness: different 
channel types exhibit a wide range of response to voltage 
while remaining functional under many mutations that 
cause their response to change, and yet they all use a 
common architecture. On these grounds, useful insights 
are expected from analyzing a mesoscale physical model. 



A. Energetics of 'idling' voltage sensors 

The mobile charges of the VS model lie within the 
electric field of the charges on the electrodes, the sta- 
tionary countercharges, and the charges induced in the 
dielectrics. We model those mobile charges as parts of 
a solid body with two degrees of freedom of solid-body 
motion: translation along the S4 axis and rotation about 
that axis. Energy is computed on a grid over this config- 
uration space [Eq. PT|) ]. From the electrostatic energy 
map, a partition function is constructed [Eq. dlSp ]. and 
from the partition function, statistics of random variables 
[Eq. P7)) ] are predicted. The energy map thus defines 
the expectations of observables such as the gating charge 
displacement (Fig. [5]) . This energy map is computed for 
each tested value of applied voltage. These maps and 
the resulting partition functions are outputs of the model 
system. 

Pseudo-color maps of electrostatic energy in the two 
dimensions computed for applied voltages of mV and 
— 100 mV are shown in Fig. [S] The translational and 
rotational origins correspond to the central S4 position 
shown in Fig. [2] for the a helix (change in energy is rel- 
ative to those origins, where Q is zero). Because of the 
symmetry of these models, the map for -|-100 mV (not 
shown) is a mirror image of the map for —100 mV. 

For a membrane voltage of mV, the energy maps of 



both helical models reveal a trough bounded on all sides 
by regions of substantially higher energy. The energy 
trough runs in the direction of proportional translation 
and rotation for the a-helical S4 model, but it runs in 
the direction of simple translation for the 3io-helical S4 
model. The trough in the map of each model follows 
the countercharge arrangement — an arrangement cho- 
sen for each model to allow periodic interactions during 
S4 motion of the S4 charges with the countercharges. In 
the a-helical model, the S4 charges and countercharges 
are aligned on parallel super-helices, whereas in the 3io- 
helical model the S4 charges and countercharges follow 
straight lines parallel to the helical axis. The lowest elec- 
trostatic energies of the a-helical S4 segment trace the 
path of a screw, whereas the energies of the 3io-helical 
S4 segment trace the path of a piston. 

The regions outside the energy trough for the 3io he- 
lix have energies about three times as large as those of 
the a helix. The electrostatic confinement of the 3io he- 
lix is stronger than that of the a helix. The strength of 
confinement correlates inversely with the separation of 
charges in the two geometries. The cluster of S4 charges 
and countercharges is more spread out in space in the 
a-helical than in the 3io-helical geometry due to the an- 
gular separations of charges in the a-helical S4 segment. 

The electrostatic energy trough tends to anchor the 
sliding-helix in a transmembrane configuration. The S4 
charges that, for a given configuration, dwell in the re- 
gion of small polarizability are balanced in this model by 
countercharges located in that region. This balance is 
maintained over the range of S4 travel where equivalent 
amounts of S4 charge and countercharge overlap in the 
region of weak dielectric (see the Supplemental Material, 
animation 1 [2y| & Fig. IH] below). A second essential 
element of balance concerns the transit of S4 charges be- 
tween the less polarizable gating canal region and the 
more polarizable vestibule and bath regions. Any energy 
change associated with the transit of an S4 charge on one 
side is approximately balanced on the other side by the 
opposite transit of an S4 charge. On the other hand, the 
energy trough generated by the electrostatics of the mod- 
els is too shallow by itself to ensure long-term stability of 
the S4 configuration. Interactions beyond those included 
in the model (such as linkages to adjacent transmem- 
brane segments or hydrophobic effect of the uncharged 
S4 residues) are necessary for long-term stability. 

To inspect the energetics more closely, we construct a 
one-dimensional energy profile for S4 translation by com- 
puting for each translational position the expectation of 
the electrostatic energy in the rotational degree of free- 
dom using the rotational partition function [Eq. (I18[) ]. 
We refer to this kind of energy profile as a 'transla- 
tional energy profile' for short. FigureslHl^E) &|ni^F) show 
the translational energy profiles for two applied voltages: 
mV {solid lines) and —100 mV {dashed lines). The 
profiles at mV are quite uniform over the translational 
extent of the energy trough (a consequence of the cho- 
sen spacing of countercharges). At a membrane voltage 
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FIG. 6. Energetics of voltage sensing. (A, C, E) a-helical S4 segment; (B, D, F) 3io-helical S4 segment. Pseudo-color maps: 
electrostatic energy for applied voltages mV (A, B) and —100 mV (C, D). (E, F): Expectation of energy over the rotational 
degree of freedom versus translation. Applied voltage Vm = rnV (solid line); Vm = —100 mV (dashed line). 



of —100 mV, the energy profiles are tilted in favor of 
more intracellular positions. A well-defined energy min- 
imum is found at a position about —1 nm inward from 
the central position of the a-helical S4 segment, whereas 
a broad minimum spread between —0.7 and —1.2 nm is 



found with the 3io-helical S4 segment. 

The energy profiles in Fig.|6]do not resemble the profile 
of an ion embedded in a lipid membrane — the latter pro- 
file has a high barrier in the center of the weak dielectric 
|28l [201 . Instead, in these models the charged section 
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FIG. 7. Energetic consequences of the countercharges. Trans- 
lational energy profile of tfie Q-lielical model of Fig. [5j A) with 
countercharges present (dashed line) or deleted (solid line); 
Cp = 4. The dotted curve is computed with countercharges 



deleted and e„ 



16. 



of the S4 helix can travel with an almost level energy 
over the translation range where S4 charges overlap with 
stationary countercharges. To examine the contribution 
of the countercharges to this result, we re-compute the 
energy profile for the a-helical S4 segment with the coun- 
tercharges deleted from the model (Fig. [7|). Deletion of 
the countercharges converts the energy trough seen with 
the standard model [dashed line) into a broad barrier 
(solid line). That barrier is reduced but not inverted 
by increasing the VS dielectric coefficient Cp from 4 to 
16 (dotted line). For comparison, a sliding-helix model 
without countercharges using gp = 10 has been analyzed 



previously by Grabe et al. 13 1 



Over the range of translation where the countercharges 
produce an energy trough [dashed line), the deletion of 
the countercharges [solid line) has a rather small effect 
(note that we plot energy relative to the central position) . 
Movement of the S4 helix encounters little energy varia- 
tion in this region of translation because the amount of 
S4 charge present in the domain of weak dielectric does 
not vary: as S4 charge enters on one side, S4 charge 
leaves on the other side. The energetics are not favorable 
for these S4 positions (as the charges are not balanced), 
but they are rather uniformly unfavorable until S4 trans- 
lation exceeds ^1 nm from its center position. If the 
travel of the S4 helix is restricted to ±1 nm in a bio- 
logical channel by means other than the countercharges 
deleted here, then even an S4 segment without balancing 
countercharges could perform as a voltage sensor. 

The variations of energy are small over the traveled 
range of translation in Figs.lHl^E) &|ni[F). The restriction 
in the model that the S4 domain and its charges must 
move as a single solid body might be expected to 'syn- 
chronize' periodic interactions among charges and coun- 
tercharges, leading the energetics to express several dis- 



tinct barriers and wells. The chosen spacing of the coun- 
tercharges, however, is enough to prevent the emergence 
of such a pattern. Additional degrees of freedom are 
thus not a prerequisite for smooth S4 travel (examples 
of such degrees of freedom are the possible flexibility of 
the individual charge-bearing S4 residues or changes in 
configuration of the helix between the a and 3io forms, 
[M2O). 

Although the energy profiles of the two helix forms 
are similar, the small differences between them are suf- 
ficient to produce substantial differences in the relations 
between displaced charge and voltage (Fig. [5]). Because 
of the small size of these energy differences, contributions 
to the Hamiltonian not included in the model (in partic- 
ular contributions arising by coupling of the VS to other 
parts of the channel) could override the differences be- 
tween the S4 helix versions seen in Fig. [S] and Figs.[BjE) 

&[S1;f). 

Additional axes of variation for these models with con- 
sequences for voltage gating include the geometry of the 
gating pore region, the dielectric coefficient of the protein 
region, the distribution of countercharges, and the effects 
of additional surface charges. We have partially explored 
these spaces [30] . The robustness of the mechanism ex- 
plored here in the face of pathology due to mutation can 
also be explored with this approach, given the computa- 
tional tractability of these mesoscale models. 



B. Energetics of '^vorking' voltage sensors 

The essence of a voltage sensor is that it can do ex- 
ternal work when the membrane voltage is varied. This 
work can then be applied, for instance, to reconfigur- 
ing the channel between conducting and non-conducting 
configurations (gating). In order to determine the work 
that a voltage sensor model might produce, we simulate 
models with a translation-dependent external workload 
included in the Hamiltonian. The potential energy field 
in which the S4 domain translates and/or rotates then 
comprises the electrostatic potential energy described by 
Eq. ((TT|) plus the potential energy Wl due to the load: 



W = Wl+W2 + WL■ 



[19) 



For our examples, we use a hypothetical load field pro- 
ducing a constant force opposing the inward translation 
of the S4 segment and hence a load potential that varies 
linearly with translation [Fig. [SjA)]. This describes the 
energetics if, for instance, the gate of a hypothetical chan- 
nel with a single VS resists closing (by inward movement 
of the S4 helix) with a constant force. 

The presence of this load alters the relation between 
the mean displaced charge and voltage by shifting the 
curve to negative voltages [Fig. [DJB), compare Fig. [5]. 
In the case of the a-helical S4 model, the charge-voltage 
curve becomes quite similar to that observed in the ex- 
periment with full channels [symbols in Fig. [SfA)]. 
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FIG. 8. Simulations of voltage sensors under a hypothetical load. (A) Component of the Hamiltonian representing the load [see 
Eq. (|19[l]. Expectations of displaced charge (B), the electrostatic energy due to the applied voltage (C), and the energy due to 
the interaction with the load (D). Solid lines: a- helical model; dashed lines: 3io-helical model. AH expectations of change in 
work are relative to the work for the configuration {z, <j)) = (0, 0) (which is identical for all V-m, see Sec. lIIIB| . 



In order to assess the relationship between the elec- 
trical work picked up by the VS and the external work 
delivered to the load, we consider the expectations of 
AM^2 [Fig.HC)] and AW^l [Fig.ISjD)]. The expectation 
for /SW2 has a simple relationship to the expectation of 
displaced charge: 



(AT4^2) = -^Q 



I J ^ra-^ij "~" y^l •'"'^ni 



(20) 



The work exchanged between the VS and the load 
((AWl)) varies most strongly with voltage where the 
variation of displaced charge is maximal, as expected if 
the load depends on S4 position. Two points on the axis 
of applied voltage are associated with polarity changes 
of displaced charge and/or of components of the Hamil- 
tonian in Fig. [5] (not marked): Vi for (Q) (B), (AVK2) 
(C), and (AVFi) (D); and I/2 = for (AW2) (C). In the 
voltage region V^a < Vi, the applied voltage drives the 



S4 helix inward while external work is done against the 
load. In the region Vm > V2, both the applied voltage 
and the work done on the S4 helix by the load drive 
the S4 helix outward ((Q) > 0). In the intermediate re- 
gion Vi < V^n < V2, the S4 helix is also driven outward 
{{Q) > 0). Here the applied- voltage and load compo- 
nents of the Hamiltonian have opposite polarities: out- 
ward S4 motion prevails because the load prevails over 
the opposing effect of applied voltage. 

The translation range traveled by the S4 helix is lim- 
ited in these simulations by the electrostatics of the VS 
rather than by the hypothesized load. The limits are set 
by the electrostatic self-energy component of the Hamil- 
tonian, AWi (Fig.m Kn = 0). Therefore, the work that 
the VS can do on the load at large (positive or negative) 
voltages saturates [Fig.ISfD)], whereas the (negative) en- 
ergy contributed to the VS by the applied field continues 
to increase in magnitude [Fig.[51[C)]. Biological K+ chan- 
nels open with very low probability (< 10~^) at large 
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negative voltages with no indication of a saturating min- 
imal probability |8l-[lG|. This may indicate that the load 
imposed on the VS domains by the 'gate' of the channel 
actually limits S4 travel at large negative voltages. In 
contrast, at large positive voltages the open probability 
of the channels saturates at levels well below 1, suggest- 
ing a saturating amount of work that the VS domains 
can do on the gate. 

Figure |8] shows the simulation results for both the a 
and 3io configurations of the S4 helix. The differences in 
voltage responsiveness observed in the simulations with 
the idling sensors are also found under the hypothetical 
load that we test. This observation, however, cannot be 
generalized. The differences between these two forms of 
VS result from the associated electrostatic self-energies 
AWi (as seen in Fig. ^. In the Hamiltonian of the sys- 
tem under load [Eq. ^T9\i ]. characteristics are determined 
by the sum of the self-energy term Wi and the load term 
Wl- Hence, it is in principle possible that the load trans- 
forms the characteristics of voltage sensing observed here 
for the 3io-helical VS model into characteristics indistin- 
guishable from those of the a-helical VS model or vice- 
versa, or a load may even result in behaviors different 
from those seen in either simulation. Experimental ob- 
servations on displaced charge therefore cannot be inter- 
preted in terms of models that include only a voltage 
sensor. These observations, like those on gating, require 
a model of the full channel to be analyzed. 

The mesoscale model system that we have presented 
for a single VS is extendable using the approach described 
in this section. The Hamiltonian of the system can be 
extended by terms describing the energetics of a chan- 
nel comprising four voltage sensors, a gating domain, 
and their coupling. The number of degrees of freedom 
(e.g., two for each voltage sensor) remains computation- 
ally manageable, so that a mesoscale model of the full 
channel can be simulated in order to obtain insight into 
the cooperation of its parts. 

We have analyzed here some equilibrium properties of 
a VS domain in which the S4 helix is relocated as a solid 
body. A refined model augmented with dynamics will 
likely have to include a second scale describing motions 
of side chains, possible deformations of the S4 helix, in- 
teractions with the bath solutions, and associated dissi- 
pative aspects. General energetic variational approaches 
developed for hydrodynamic systems of complex fiuids 
31| may provide a multi-scale approach to channel dy- 



Holohean, Dr. Peter Larsson, and Dr. Karl Magleby for 
helpful discussions. 

Appendix: Charge distributions of VS models 

The sliding helix in our VS models is a microscopic 
voltage sensor. Therefore, we describe the VS by its sta- 
tistical ensemble behavior. In this appendix, stochastic 
VS behavior is visualized for several models presented 
in the paper, both as figures for a fixed applied voltage 
(—100 mV; Figs.|n]&[ini) and as associated Supplemental 
Materials, animations 1 ^26,] & 2 [23] with voltage increas- 
ing uniformly over time from —100 mV to -1-100 mV. 

The figures show two stochastic aspects of VS behav- 
ior: (1) the mean positions of the S4 charges [marked by 
blue (dark gray) balls], and (2) the charge density distri- 
bution of S4 charge [represented by a blue (gray) cloud 
with a color intensity proportional to the charge density 
there] . A high density of color marks the locations where 
the S4 charges dwell frequently, as opposed to their mean 
positions. 

The expectation of position for each charge is com- 
puted from Eq. (IT7)l using the positions r^ for the charge 
Qk as the random variable X, over the partition func- 
tion with translational and rotational degrees of freedom. 
Since the helix behaves as a solid body, the helix position 
r fixes the positions r^ for the charges qk ■ That relation- 
ship allows us to define the partition function, the energy 
function, the configuration probability and our measures 
for the positions r^ in terms of the respective function 
for the helix position r. For the expectation of position 
for each charge: 



(r.) 






2. J 



^ijk-Tij 






Wij/keT 



(A.l) 



where Wij is the work needed to construct configuration 
ij, Pij is the probability of that configuration, and Vijk 
is the position of charge k in configuration ij. 

Likewise, the distribution of charge can be computed 
by applying the partition and energy functions in terms 
of the positions rfc of charges qk- The charge density 
z(r) is then the sum over all charges of the probability of 
each charge being located at r, multiplied by its valency 
in units of Co, and normalized: 
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where 6{r,rijk) is the discretized delta function (1 if we 
are treating Vijk as the same location as r for visualization 
purposes; otherwise 0). 

The color representations for the animations are pro- 
portional to z(r), normalized to the highest charge den- 
sity at that frame's potential. 
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FIG. 9. Standard a-helical model: Position and distribution of charges for model a (1) in Table [T]^C) with ep = 4. Red (light 
gray) symbols represent fixed negative countercharges, blue (dark gray) symbols represent the mean position of —1/3 eo on S4 
arginines, and blue (gray) shading represent the relative probability of ne gati ve charge at a given position. This is the model 
used in Fig. [5jA) . See the associated Supplementary Material, animation 1 2©] for the behavior of this model of the VS domain 
over the range of transmembrane potentials from —100 mV to +100 mV. 




FIG. 10. Standard 3io -helical model: Position and distribution of charges for model 3io in Table [T](C) with ep = 4. See 
the description of Fig. [9] for further details. This is the model used in Fig. [5jB). See the associated Supplemental Material, 
animation 2 [23] for the behavior of this model of the VS domain over the range of transmembrane potentials from —100 mV 
to +100 mV. 
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